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Abstract 



This paper investigates Monte Carlo techniques for construction of compact 
wavefunctions for the internal atomic motion of the D3O''' ion. The polar- 
ization force field models of Stillinger, et al. and of Ojamae, et al. were 
used. Initial pair product wavefunctions were obtained from the asymptotic 
high temperature many-body density matrix after contraction to atom pairs 
using Metropolis Monte Carlo. Subsequent characterization shows these pair 
product wavefunctions to be well optimized for atom pair correlations de- 
spite that fact that the predicted zero point energies are too high. The pair 
product wavefunctions are suitable to use within variational Monte Carlo, in- 
cluding excited states, and density matrix Monte Carlo calculations. Together 
with the pair product wavefunctions, the traditional variational theorem per- 
mits identification of wavefunction features with significant potential for fur- 
ther optimization. The most important explicit correlation variable found for 
the DaO"*" ion was the vector triple product roDi-{'<^OD2^^0D3)- Variational 
Monte Carlo with 9 of such explicitly correlated functions yielded a ground 
state wavefunction with an error of 5-6% in the zero point energy. 
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INTRODUCTION 



Flexible and dissociative models for simulation of liquid water have often been used 
in classical statistical mechanical studies of aqueous materials. They present technical 
advantages for carrying-out simulations and are of-the-essence where dissociation of wa- 
ter molecules is necessary to the chemistry being studied. Studies of clusters suggest that 
quantum mechanics plays a non-negligible role in proton transfer .0 However, directly incor- 
porating quantum mechanics via discretized path integral approaches may require orders of 
magnitude larger computational effort than the classical problem. This Report investigates 
constructing simple intramolecular wavefunctions for aqueous species. These wavefunctions 
might be used in Monte Carlo simulations of aqueous solutions with the same computational 
techniques used to treat flexible simulation models classically. More specifically, our goal is 
to determine the complication and accuracy to be expected in constructing wavefunctions 
that might be transferrable and useful in related studies of aqueous solution chemistry. 

We take the deuterated hydronium ion, D3O+, as a specific example species. We choose 
this ion because of the substantial current interest in proton exchange processes in water and 
in the spectroscopy of charged clusters of water molecules;lllil&i§ we choose the deuterated 
case to avoid the interesting complications of spin restrictions on the wavefunctions. In 
addition we will study here only states of zero angular momentum, J=0. Our goal will be 
to obtain a simple intramolecular wavefunction for this species in a reasonably organized 
fashion. 

A traditional exhaustive expansion of the wavefunction into a basis will produce satis- 
factory vibrational energies with effort but does notproduce a simple result for the wave- 
function. Less traditionally, diffusion Monte Carloc3 will produce satisfactory vibrational 
energies for low energy states with effort but not a simple wavefunction for other uses. 

In contrast, the historical work of McMillancil, a variational Monte Carlo calculation for 
the ground state of liquid He'', does directly center on the construction of a simple wavefunc- 
tion. The Monte Carlo character of this technique is limited to the primitive but essential 
task of evaluating integrals with a many-body wavefunction. This approach permits sim- 
ple but sophisticated descriptions of correlation and can include a limited number of excited 
states. The possibility of simultaneous treatment of excited states provides an additional av- 
enue for description of correlation. In fact, the density matrix Monte Carlo approach0&0 
of recent years can be regarded as a generalization of the McMillan calculation; the (3=0 
circumstance is precisely the result that we view as the McMillan calculation with excited 
states. 

The full power of the density matrix Monte Carlo approach can be used with rough 
initial estimates of the wavefunctions sought. The Monte Carlo calculation systematically 
improves the computed energy levels. However, since our goal is to derive a wavefunction 
useful in other contexts, we emphasize finding and utilizing the best simple function we can 
before the density matrix Monte Carlo procedure takes over. 

We will initially use pair-product wavefuctions for atomic motion. The components of the 
pair-product are obtained by diagonalizing an approximate density matrix. We emphasize 
that initial model vibrational wavefunctions need not be orthonormal. 

It is found (below) that these initial pair-product wavefunctions capture two body cor- 
relations fairly well. Using them as a basis in McMillan and density matrix Monte Carlo 
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calculations incorporates further correlation effects. To obtain more compact wavef unctions, 
the variational principle is used to identify the important many-body correlations. These 
effects are included into the ground state wavefunction and yield a significantly improved 
simple wavefunction for the ground state. A few of these correlated wavefunctions can be 
used as the basis in a McMillan calculation yielding reasonably accurate wavefunctions. 

MODEL WAVEFUNCTIONS FOR ATOMIC MOTION 

We develop the vibrational wavefunctions analogously to development of wavefunctions 
for liquid He. We initially seek a pair product form 

•^0 = n (PoDiroD) n (pDoirDD) ■ (1) 

OD pairs DD pairs 

In this case, we expect the functions (poD to be local mode orbitals and the functions (pco 
to serve as correlation functions. 

The casual interpretation of these functions as orbitals will serve to distinguish vibra- 
tional excitations. The method sketched here is based upon the observation that the eigen- 
functions of the asymptotic /3 — > density matrix are suitable initial estimates of the pair 
functions sought. This method has the additional advantage that harmonic analysis of a 
potentially rough energy landscape is not required. 

We start with a function of this simple form for several reasons. First, a similar distri- 
bution is implied when flexible and dissociative pair simulation models of water are treated 
classically. Second, elements of such a form might be approximately transferrable to OD and 
DD joint distributions in other settings - this would be the natural assumption for simulation 
calculations. Third, such a form would be similar, though not the same as a ground state 
vibrational wavefunction for a harmonic system. The distinctions are that normal modes 
coordinates are not used, and that functions employed can attempt to treat anharmonic sys- 
tems by deviating from Gaussian form where the potential is anharmonic. The remainder 
of this paper characterizes such wavefunctions and investigates some improvements. 

Natural orbitals for general oscillators in one dimension 

The natural orbitals of a density matrix {x'; l3\x) =< x'\e~^^\x > are introduced by 

J {x'; f3\x)M^')dx' = e-^^'-M^) ■ (2) 
In fact the asymptotic (3 —>■ limit of the density matrix is simpleS 

{x'; (3\x) = {m/2T^ph^f/^ exp{-m(a; - x'f/2ph^ - p{y{x) + V{x'))/2} . (3) 

Sethia et al. showed that though the density matrix used is approximate, satisfactory orbital 
functions are obtained.^ This form for the thermal density matrix is more than merely an 
approximation. For small /5, it is asymptotically correct. Therefore, it is satisfactory if the 
application permits a small (3. The subsequent developments exploit this point. 
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Contraction to Pairs for the many-body Thermal Density Matrix 



For higher dimensional problems, we can exploit the approximate density matrix to 
produce the necessary reasonable pair functions after tracing-out all other degrees of freedom. 
This tracing-out will be generally possible on the basis of classical Monte Carlo techniques. 

Consider a system composed of N particles located at (ri, r2 . . . ftv)- The N-body high- 
temperature thermal density matrix is 

N 

(n, r2 . . . r^; P\r\, r'2 . . . rV) = H imk/27Tl3hY^ exp{-mfc(rfc - r'kf/2(3h^} 

k=l 

X exp{-/3(l^(ri, . . . r^) + \/(r'i, r'2 . . . rV))/2} . (4) 

We now focus on the {ij} pair of particles. For all other particles, we restrict the density 
matrix to the diagonal. For example, taking {ij}={12} we obtain 

(ri, r2 . . . r^v; f3\r'i, r'2 . . . r^v) oc H (^fc/Svr/J/i')'/' exp{-mfc(rfc - r'k)y2(3h'} 

k=l,2 

X exp{-(3{V{r,, r2 . . . r^) + ^(r'l, r'2 . . . r^))/2}. (5) 

Now for the {12} pair we transform to the {12} center of mass and relative coordinates: 
R = (miri + m2r2) / (mi -|- 1712) and r = ri — r2. 

(R + m2r/M, R - mir/M . . . r^v; /3|R' + mgr'/M, R' - mir'/M . . . r^v) oc 

exp{-(3{V{K + mar/M, R - mir/M . . . r^v) + ^^(R' + mar'/M, R' - niir'/M . . . rjv))/2} 
X {M/27r(3n''f/^ exp{-M(R - R') 72/3^' }(/x/27r/?;i2) 3/2 exp{-/i(r - r')V2/?^'}, (6) 

where M = mi + m2 and /i = 71117712/ M. 

Note that R is not the molecular center of mass and that the potential generally does 
depend on R. However, in the interest of simplicity and in view of the form sought Eq. 
(0), we will trace-out the R dependence also. Our final step will be to use classical Monte 
Carlo techniques for the integrations required to bring Eq. into the form of Eq.(^ 
for the relative coordinate r. Thus we sample configurations with the probability density 
exp{—pV{ri, r2 . . . rjv)} and we estimate the kernel 

where 

(5V(R + mar'/M, R - ttht'/M . . . r^) = 
\^([miri + 7n2Y2 + m2r'] /M, [miri + m2r2 - mir'] /M . . . r^) - V{ri, r2 . . . r^) . (8) 

K{r, r') is symmetric. For the case considered here, it does not depend on the angles that r 
or r' make with the laboratory fixed coordinate system. Thus, this kernel depends only on 
r, r', and an angle. Since we here focus on radial functions, we only need the kernel after 
having averaged over that polar angle. 

Note that for determination of the orbitals only, the normalization of the kernel K{r, r') 
is not significant. To within an unimportant normalization constant we can evaluate the 
required kernel through the following procedure: 
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1. Draw configurations (ri, r2 . . . Tat) from tlie probability distribution proportional to 
exp{—pV{ri, r2 . . . rjv)} utilizing the Metropolis Monte Carlo algorithmiJS. 

2. For each configuration, choose an r' from a grid and its corresponding spherical angles 
6' and (p' using quasi-random number series. Weight each configuration (ri, r2 . . . ttv) 
and the corresponding r' by 

3. Perform a final integration over angles 

K{r,r') oc f d9 ^ d9' d<p' ^"^0 sin^ K{v,v') sin^' (10) 
Jo Jo Jo Jo 

6,(j),6', and 0' are the spherical, angles corresponding to r and r'. For this case of 
J = 0, it would have been sufficient to integrate over the angle between r and r'. 

With this kernel in hand, we solve the one dimensional equation 

/ K{R,R')R'^(f){R')dr' = X^{R) . (11) 



A is not e~^^ but is proportional to it. 

This approach is not limited to the circumstances that one member of that pair is a 
massive molecular center. This can be applied also for DD pairs in the DsO"*" molecule. 
However, for the D3O"'" molecule in particular, it is natural to regard the OD functions as 
local mode orbitals and the DD functions as providing a subsequent account of correlations. 



D3O+ PAIR FUNCTIONS 

The functions (poD and (pDD are found for the second version of the Ojamae-Shavitt- 
bmger (OSS) potential! using (3 ^=0.01 Hartree and a sample size of 500,000 (ri, r2 . . . r^r) 
configurations. This energy parameter is 6.3 kcal/mole (more than 3000 K), a value above 
the inversion barrier of 4.4 kcal/mole. Smaller values for jS are problematical because of the 
fragmentation of the DsO"*" molecule at high temperatures. 

Figure |l| shows (puu obtained. The additional approximation of fixing one of the 
deuterons in the reference bond has the significant effect of narrowing the pair function. 
As expected, the effect is not as dramatic for fixing an oxygen atom. The co-linear assump- 
tion of fixing the orientation of the reference bond is not significant. 

These functions were also obtained for the Stillinger, Stillinger, Hogdgon (SSH) 
potential.0 A comparison of features of SSH and OSS is presented in Table |. The higher 
(more realistic) binding energy for the additional proton in the SSH case and the higher 
inversion barrier of this potential permits a smaller f3. In that case, the pair functions are 
found for two values of namely 0.01 Hartree and 0.01677. . . Hartree. The difference 
in the "frequencies" associated with the ground and first excited state OD "local mode" 
functions changed from 3150 cm~^ to 1987 cm~^. The corresponding differences in DD fre- 
quencies changed from 1243 cm~^ to 837 cm~^. In spite of these large frequency changes, 
the pair functions, seen in Figure ^, are qualitatively similar. 
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ACCURACY OF THE PAIR PRODUCT WAVEFUNCTIONS 



The numerical arrays describing the OD and DD functions were fit to cubic sphnes with 
first derivatives set to zero at the end pointsU and the pair product wavefunctions of Eq. (|l]) 
were obtained. Several steps were taken to assess the quality of these wavefunctions. For 
the OSS and SSH potentials, the variational Monte Carlo0 zero point energies of the pair 
product wavefunctions are 0.0321(2) and 0.0367(6) Hartrees, respectively. Simple diffusion 
Monte CarloS zero-point energies for the potentials are 0.02516(3) and 0.03042(3) Hartrees. 
(All errors are two standard deviations unless otherwise noted.) The energy for the pair 
product function is between 20% and 30% too large. This function was further optimized 
by introducing parameters to move (m) and scale (s) the pair functions; (pivnew) = <P{{''^oid + 
m) I s) . No significant lowering of the energy was found. 

The quality of optimization of these pair functions can be analyzed directly. The usual 
variation principle can be expressed as 

((51n^)Ei(l..iV))|^p - ((51n^))|^|2 (Ei(l..iV))|^|2 = , (12) 

where Ei^ = "if'^H"^ is the usual local energy function and the subscripted brackets (. . .)|^|2 

indicate an average with probability density |\E'| . For a pair product wavefunction this form 
implies 



P?^^W)|^|2(^L(l..iV))|^|: 



, (13) 



where p^J{r) is the 'ju pair joint density operator. 

This stationary requirement is satisfied for the exact ground state wavefunction since 
the local energy is spatially constant. Deviation of the left side of Eqn. |I3| from zero will be 
defined as optimization deficit. 

Figure ^ shows the optimization deficit for the OSS DD pair functions and the DD VMC 
radial distribution function. The OD functions show similar low optimization deficits; these 
pair functions lead to accurately uncorrelated pair densities and local energies and are thus 
well optimized. Note that the largest deficits are at the tails of the distribution. Thus further 
optimization would require significant effort because of the difficulty of getting statistically 
significant data for those infrequent events. 



CORRELATION BEYOND PAIRS: VARIATIONAL AND DENSITY MATRIX 

MONTE CARLO 

We can use model excited state functions to improve the ground state and get a few 
low lying states at the same time. A concise expression of these approaches is obtained by 
focusing on the density matrix. Suppose we have a number of approximate wavefunctions 
|$j),z = 0, . . . , A^. By Monte Carlo procedures detailed elsewhere^'El, we estimate the 
matrix elements 

p,,.(/3)=<<l>,|exp{-/?iJ}|$, >, (14) 
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and 



H,0) =< ^,\Hexp{-/3H}\^^ >, (15) 

for P not too large. As P gets large, information is lost for highly excited states, so this 
approach depends on location of an intermediate range of f3. H and p are then diagonalized 
simultaneously. The need not be orthonormal initially since Monte Carlo methods 
will be used to perform the required integrations. The /3 — limit may be viewed as a 
direct extension of the McMillanEH approach. The trace of p^^ provides a natural importance 
function for the Monte Carlo calculation. 

To apply these methods to our problem, an ordered set of excited state functions is 
constructed. The ratio of eigenvalues of the n^^ and 0*^ eigenfunctions in Eq. |Tl| are used 
to assign relative energies to the pair product functions. Model excited state functions are 
built-up by linearly combining functions representing higher levels of excitation. 

These functions were used in density matrix Monte Carlo calculations. In its /3 = 
limit, density matrix Monte Carlo mixes the lowest energy functions much like configuration 
interaction in electronic structure calculations. Tables H and |T| show the energies obtained 
in these variational calculations for OSS and SSH potentials, respectively. The guiding 
function used was \l/triai(Q') = {j2i-Lo ^ii^)"^^ with M=5 and = \. z/ of | is used instead 
of I so that higher states than M are efficiently sampled. Alternatively, M could be made 
larger. However, the cost of making M larger is higher than the cost of making = \- 
Independent estimates for the zero-point energy of 0.02516(3) Hartree (5522(7) /cm) and 
0.03042(3) Hartree (6676(7) /cm) are found using diffusion Monte Carlo for the OSS and SSH 
potentials, respectively. Inclusion of a small number of excited states substantially improves 
the predicted zero-point energy; the 20-30% error with the pair product wavefunction is 
reduced to about a 10% error with inclusion of nine states. However, the improvement of 
this prediction with increasing numbers of excited states is slow. As shown in Tables H 
and m, inclusion of 94 and 99 states still results in a 4-6% error. 

The variational energies are significantly improved by density matrix Monte Carlo as 
f5 increases from zero with 9 excited state functions. Including more functions accelerates 
the convergence, (see Figure ^) Tables |IV| and |V| show the energy estimates for lowest 5 
symmetric states for OSS and SSH respectively. The average energies and their standard 
deviations as a function of P were used in a Marquardt's fit to the function A + Be~^^^. 
Error bars were obtained by dividing the configurations used to calculate the matrices into 
smaller sets, diagonalizing these sets, obtaining the standard deviation for their eigenvalues, 
and dividing it by the square root of the number of sets. Assuming that the distribution of 
errors is Gaussian, the ellipsoid representing 90% error can be plotted.0 The errors reported 
in Tables |I^ and are 90% confidence intervals. Some data was not well approximated 
by a decaying exponential. In this case, the point in the plateau with the smallest standard 
deviation is given. For this data which is designated by a *, 2 standard deviations are given 
in parenthesis. 

Because it is relevant to assessment of potential energy surface models, we note in passing 
that for the OSS potential the excited vibrational state with a node roz3i-(roD2XroD3)=0 
has an energy of 0.02522(4) Hartree as determined by a fixed node diffusion Monte Carlo 
calculation. This is about about 15(24) /cm higher than the zero point energy. A model fit to 
experiment by Sears et al. yielded a gap of 15/cm. However, the OSS potential overestimates 



7 



the inversion barrier and, therefore, we expect it to underestimate this sphtting. Beyond 
this comparison, our results for both the OSS and SSH potential differ signicantly from the 
results of Sears et al. For example, first totally symmetric excited states are 794(10)/cm 
and 731(ll)/cm above the ground state for OSS and SSH respectively. The first totally 
symmetric excited state in Sears' fit was 453/cm above the ground state.S 

EXPRESSING CORRELATIONS BEYOND PAIRS MORE COMPACTLY 

The variational principle Eq. |l^ and the optimization deficit can be used to search for 
structural variables that are correlated with the local energy. Such variables are principal 
candidates for construction of explicitly correlated wavefunctions. 

Consideration of several natural possibilities for such structural variables identified the 
coordinate Vi=v o or 002^^003) as more correlated with the local energy of the pair prod- 
uct wavefunction than any other combination considered. The distribution of u and the 
optimization deficit (see Fig. ^ with this pair product wavefunction suggested the form 

^n-body = e-'^(«-^)Ve-'^("+^)' . (16) 

The pair product wavefunction with moving and scaling parameters times this new mul- 
tiplicative many-body factor was optimized using Variational Monte Carlo. Zero point 
energies of 0.02838(4) Hartree (6229(9) /cm) and 0.03326(2) Hartree (7300(4) /cm) are ob- 
tained for the correlated wavefunctions with the OSS and SSH potentials, respectively. This 
is nearly as good as the VMC result with mixing of 9 states. The optimized parameters are 
shown in Table 0. The new distribution and optimization deficit for the OSS potential are 
shown in Fig. ^. 

Using pair product wavefunctions with this correlation piece as the basis in a McMillan 
VMC calculation with 9 functions reduced the ground state energy to 5860(2) /cm and 
7041(3) /cm for OSS and SSH potentials, respectively. These energies contain 5-6% error 
with respect to the DMC ground state energies. 

CONCLUSION 

A pair product wavefunction for D3O+ has been obtained from an approximate density 
matrix. This pair product wavefunction captures two body interactions fairly well. Varia- 
tional Monte Carlo includes correlation by mixing in excited states. However, it takes about 
99 states to get to an error of 4-6%. Density Matrix Monte Carlo gets accurate energies 
with 9 wavefunctions. However, it does not provide a simple wavefunction. 

A compact ground state function is obtained by using the variational principle to iden- 
tify the most significant many-body terms and including it directly into the wavefunc- 
tion. For the DaO^ ion the most significant such variable was the vector triple product 
roDi-(roD2XroD3)- Variational Monte Carlo with 9 of these correlated functions yields a 
ground state wavefunction with an error of 5-6% in the zero point energy. 
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FIG. 1. The <?!>dd('') obtained by fixing one of the deuterons (hne with open circles) is sig- 
nificantly narrower than that obtained from unconstrained sampling (line with diamonds). The 
co-linear assumption of fixing the orientation of the reference bond (line with solid circles) is not 
significant. The points were fit to cubic splines. Displacements are in A . 
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FIG. 2. The (Pbd{i") obtained using l//3=0. 01677. . . Hartree (open circles on dashed hne) is 
shghtly wider than that obtained using 1//3=0.01 Hartree (diamonds on sohd hne). The points 
were fit to cubic sphnes. Displacements are in A . 
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FIG. 3. DD optimization deficit (diamonds on solid line) and VMC radial distribution function 
(open circles on dashed line) for OSS model. 
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FIG. 4. Convergence of the lowest energy state to the ground state energy as a function of (3 
for the OSS potential. The circles and squares were obtained using 9 and 17 states, respectively. 
The lines are Marquardt's fits. 
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FIG. 5. The OSS VMC distribution of u (solid line) and the optimization deficit (dashed line) 
is shown for the pair product wavefunction 
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TABLES 

TABLE I. Comparison of SSH and OSS potentials with B3LYP hybrid density functional theory 
calculations of Martin, Hay, and Pratt.^ 



Quantity 


SSH 


OSS 


ab initio^ 


H+ Binding on H20 


0.269 Hartree 


0.1695 Hartree 


0.27367 Hartree 


H30+ Inversion Barrier 


4.78 kcal/mole 


4.42 kcal/mole 


2 kcal/mole 


zero-point (Normal Modes) 


0.0339 Hartree 


0.0281 Hartree 


0.0212 Hartree 


zero-point (DMC) 


0.03042 Hartree 


0.02516 Hartree 





^ R. L. Martin, P. J. Hay, and L. R. Pratt. J. Phys. Chem. A 102, 3565 (1998). 



TABLE II. Variational Energies for the first 5 states using 9, 17, 46, and 99 functions for the 
OSS potential. An independent estimate for the zero-point energy of 5522(7)/cm is found using 
diffusion Monte Carlo. 64 trajectories generating 10,000 configurations were used in the integrals. 



States 


9 


17 


46 


99 





6198(2) 


(3083(4) 


5903(1) 


585(3(4) 


1 


7638(4) 


7213(6) 


7001(6) 


6905(5) 


2 


9588(5) 


9106(10) 


8702(12) 


8508(10) 


3 


10045(6) 


9449(12) 


8804(8) 


8672(7) 


4 


10176(6) 


9677(8) 


9269(8) 


9105(8) 



TABLE HI. Variational Energies for the first 5 states using 9, 17, 44, and 94 excited state 
functions for the SSH potential. An independent estimate for the zero-point energy of 6676(7)/cm 
is found using diffusion Monte Carlo. 64 trajectories generating 10,000 configurations were used in 
the integrals. 



States 


9 


17 


44 


94 





7241(4) 


7112(4) 


6974(4) 


6939(4) 


1 


8692(6) 


8201(6) 


8065(6) 


7968(6) 


2 


10604(6) 


10176(10) 


9790(8) 


9538(8) 


3 


11498(14) 


10746(16) 


10043(8) 


9885(8) 


4 


11927(8) 


11386(8) 


10729(8) 


10484(6) 
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TABLE IV. Density matrix Monte Carlo energies for the first few states using 9 and 17 symmet- 
ric excited state functions for OSS. One million configurations where used in density matrix Monte 
Carlo. The time step was 1 / Hartree. Error bars are 90% confidence intervals. An independent 
estimate for the zero-point energy of 5522(7)/cm is found using diffusion Monte Carlo. 



States 


9 


17 





5525(4) 


5513(5) 


1 


6297(10) 


6307(5) 


2 


6883(22) 


6881(28) 


3 


7979(10) 


7790(38) 


4 


8083(18) 


8038(20) 



TABLE V. Density matrix Monte Carlo energies for the first few states using 9 and 17 sym- 
metric excited state functions for SSH. One million configurations where used in density matrix 
Monte Carlo. The time step was 1 / Hartree. Error bars are 90% confidence intervals. Some data 
was not well approximated by a decaying exponential. In this case, the point in the plateau with 
the smallest standard deviation is given. For this data which is designated by a *, 2 standard devi- 
ations are given in parenthesis. An independent estimate for the zero-point energy of 6676(7)/cm 
is found using diffusion Monte Carlo. 



States 


9 


17 





6661(5) 


6675(4) 


1 


7437(7) 


7406(7) 


2 


7964(54)* 


7934(78)* 


3 


8912(4) 


8918(10) 


4 


10119(25) 


9061(43)* 


TABLE VI. Parameters used for the variational ground state wavefunction. 




OSS 


SSH 


SOD 


0.911 


0.853 


SDD 


0.935 


0.895 


moD 


0.157 


0.242 




0.078 


0.100 


a 


9. 


10. 


b 


0.851 


0.911 
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